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Summary 

^^A  fpurth-order  Orr-Sommerfeld  eigenvalue  solver  for  heated 
boundary?- layers  developed  by  BMT  fctdk  has  been  implemented  on  the 
VAX-11/785  computer,  at  ARE<T«ridlngtonK-  The  solver  is  based  on 
the  compound  matrix  method  and  the  effects  of  heating  are  limited 
to  mean  viscosity  variation  across  the  boundary  layer.  Some  test 
computations  have  been  carried  out  on  flat  plates  with  uniform 
wall  temperature.  For  small  overheat,  the  results  agree  with 
those  previously  published  using  other  methods,  but  for  large 
overheat,  there  are  some  minor  differences.  The  calculation 
procedure  requires  a  large  number  of  integration  steps  across  the 
boundary  layer  to  obtain  Accurate  eigenvalues,  but  is  relatively 
insensitive  to  other  physical  and  numerical  parameters,  A  few 
computations  have  been  carried  out  using  a  full  sixth^order 
models  also  developed  by  BMP  which  takes  into  account 

fluctuation^  in  all  the  fluid  properties.  It  is  concluded  that 
the  fourth/order  method  can  be  incorporated  into  an  axl symmetric 
body  calculation  procedure.  t->  ■>*  ^  sr. _ 
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LIST  OF  SYMBOLS 

Cp  Specific  heat  capacity  (J  kg-1*"1) 

c  Dimensionless  complex  wavespeed  (*  w/aue> 

f  Dimensionless  stream  function,  defined  by  equation  (A12> 

k  Thermal  conductivity  <W  m”^K_1> 

R  Displacement  thickness  Reynolds  number  ( *  u^s'  Vv.) 

T  Absolute  temperature  (X) 

u  Velocity  component  in  x-directlon  (m  s”1) 

u  Dimensionless  velocity  (*  u/u  ) 

e  _i 

v  Velocity  component  in  y-directlon  (m  s  ) 

x  Streamwlse  co-ordinate  (a) 

y  Normal  co-ordinate  (m) 

Complex  wavenumber 

Dimensionless  complex  wavenumber  (■  »S*) 

Boundary- layer  displacement  thickness,  defined  by 
equation  (1) 

Transformed  normal  co-ordinate,  defined  by  equation  (A9) 

-2 

Dynamic  viscosity  (N  s  m  > 

Dimensionless  viscosity  («  u/M„> 

Kinematic  viscosity  (m^s-1) 

Transformed  boundary- layer  streamwlse  co-ordinate 
Density  (kg  m~3) 

Amplitude  of  disturbance  stream  function 
Angular  frequency 

Dimensionless  angular  frequency  («  w6*/ue> 

Subscripts 

e  Value  at  outer  edge  of  boundary  layer 

1  Imaginary  part  of  complex  quantity 

r  Real  part  of  complex  quantity 

w  Value  at  plate  surface 

•  Freestream  value 
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The  drag  of  an  underwater  vehicle  can  be  substantially 
reduced  by  delaying  transition  from  laminar  to  turbulent  flow  in 
the  boundary  layer.  This  results  in  higher  speed  for  a  given 
power  input  or  a  reduced  power  input  to  achieve  the  sane  speed. 
One  way  of  delaying  transition  is  by  heating  the  surface.  The 
viscosity  of  water  decreases  with  temperature,  and  this  phenom¬ 
enon  combined  with  surface  heating  results  in  a  more  stable 
boundary  layer. 

The  prediction  of  the  transition  location  on  a  heated  body 
can  be  carried  out  by  adapting  the  well-established  procedures 
for  an  unheated  body  based  on  linear  stability  theory  Cl]. 
Starting  at  the  nose  and  working  downstream,  laminar  boundary- 
layer  mean-flow  profiles  are  calculated  at  several  positions  on 
the  body  surface.  For  each  set  of  profiles,  the  spatial  growth 
rates  of  two-dimensional  disturbances  are  calculated  at  several 
frequencies  by  obtaining  eigenvalues  of  the  fourth-order  Orr- 
Sommerfeld  equation.  The  growth  rates  at  a  specific  frequency 
are  then  integrated  along  the  body  surface  and  the  location  where 

the  total  amplification  ratio  is  en  (with  n  usually  having  the 
value  9)  is  calculated.  The  most  forward  of  these  locations  is 
the  predicted  transition  location.  Although  the  method  assumes 
parallel  mean  flow  and  two-dimensional  disturbances,  it  has  been 
widely  verified  experimentally. 

The  objective  of  the  present  investigation  is  to  validate  a 
fourth-order  eigenvalue  solver  developed  by  BUT  Ltd  based  on  the 
compound  matrix  method.  The  heated  flat  plate  with  uniform  wall 
temperature  was  used  as  a  test  case,  since  the  partial  differen¬ 
tial  equations  describing  the  mean  flow  reduce  to  ordinary 
differential  equations  which  can  be  solved  very  accurately  and 
eigenvalue  computations  can  be  compared  with  previously  published 
results.  A  few  numerical  results  have  been  obtained  using  a 
sixth-order  system  of  equations  which  Includes  viscosity, 
temperature  and  density  fluctuations  and  buoyancy  effects. 


2.  THEORY 

The  equations  of  the  mean  flow  over  a  flat  plate  with 
uniform  temperature  are  derived  and  stated  in  Appendix  A  (equa¬ 
tions  (A22)  and  (A23))  along  with  the  boundary  conditions  ((A18), 
(A19)  and  (A21)).  The  solution  of  these  equations  gives  profiles 
of  dimensionless  velocity  and  temperature  and  their  first  and 
second  derivatives.  The  only  physical  input  required  is  a 
knowledge  of  the  fluid  properties  (density,  viscosity,  specific 
heat  capacity  and  thermal  conductivity),  which  are  assumed  to 
vary  only  with  temperature.  Two  sets  of  expressions  for  these 
quantities  have  been  used,  proposed  by  Lowell  &  Reshotko  C2]  and 
JCaups  £  Smith  C33.  Full  details  are  Included  in  Appendix  B. 
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The  modified  Orr-Soamerfeld  equation  requires  as  input  the 
dimensionless  velocity  profile  u/uA  and  its  second  derivative, 

®  £ 

expressed  as  functions  of  the  displacement  thlcicness  unit  y/S  , 

where  8  represents  the  usual  boundary- layer  displacement 
thickness,  qiven  by 

It  also  requires  the  dlaenslonless  viscosity  p/pe  and  Its  first 
and  second  derivatives  (which  depend  on  teaperature)  expressed  as 

functions  of  y/S*.  However,  the  velocity  and  temperature 
profiles  and  their  derivatives  output  by  the  mean  flow  calcu¬ 
lation  are  functions  of  t),  the  compressible  scaling  parameter 
(equation  (A9),  Appendix  A).  To  convert  the  mean  flow  profiles 
to  a  form  suitable  for  input  to  the  eigenvalue  calculation,  a 
numerical  procedure  and  FORTRAN- 77  computer  program  (known  as 
NMIAL514)  were  written  and  supplied  by  NMI  Ltd  (now  BMT  Ltd)  C43. 
In  a  typical  boundary-layer  calculation,  80  to  100  points  are 
normally  sufficient  to  determine  the  mean  velocity  and  tempera¬ 
ture  profiles,  whereas  the  number  of  data  points  (or  Integration 
steps)  required  for  accurate  eigenvalue  calculation  is  much 
larger  (typically  800).  Program  NMIAL514  also  carries  out  the 
required  Interpolation  to  achieve  this  using  a  cubic  polynomial 
method.  As  originally  supplied,  the  program  reads  in  profiles  at 
81  input  points  and  interpolates  them  to  801  output  points 
(equivalent  to  800  integration  steps),  but  it  has  been  modified 
so  that  both  these  numbers  can  be  varied. 

The  eigenvalue  calculation  method  solves  the  Orr-Soamerfeld 
equation 


(u  -  c>(4"  -  a?  4)  -  u"  4  +  -i-  J  p  (4""  -  2  a2  4"  +  a4  4) 

5  R  { 

+  2  p'  (♦'"  -  «2  4')  ♦  M"  (4"  +  *2  4>1  *  0  (2) 


which  has  been  modified  to  Include  viscosity  variations  in  the 
mean  flow.  The  mean  flow  is  assumed  to  be  parallel  and  the 
stream  function  representing  a  single  disturbance  is  assumed  to 
be  of  the  form 


4  (*.  7,  t)  ■  4(y>  exp  Cia(x-ct)}  (3) 

The  boundary  conditions  are 

4(0)  -  4' (0)  -  4(»)  ■  4' («■>  -  0  (4) 

A  numerical  procedure  and  computer  program  (known  as  NMIHE40S ) 
has  been  supplied  by  NMI  Ltd.  C4J  based  on  the  compound  matrix 
method  C53.  This  particular  method  was  used  to  avoid  any 
numerical  problems  caused  by  the  fact  that  the  Orr-Sommerfeld 
equation  is  stiff.  The  program  was  written  in  an  interactive 
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fora,  but  it  has  since  been  modified  to  allow  for  automatic 
calculation  of  eigenvalues.  The  eigenvalue  relation  has  the  form 

F  (a,  w,  R)  «  0  with  a  ■  ar  +  Sj  and  w  *  i3r  +  uij  (5) 

There  are  three  nodes  of  operation.  Mode  0  gives  temporal 
Instabilities.  «r  and  are  specified  (with  oi^  ■  0)  and  the 

program  calculates  u»r  and  Sj.  A  positive  value  of 
indicates  an  unstable  mode.  Mode  1  (spatial  instability)  is 
Where  wr  and  are  specified  (w^  =  0)  and  and  are 

calculated.  Mode  2  (also  spatial  instability)  is  where  and 

Wj  *  0  are  specified  and  and  ur  are  calculated.  A 

negative  value  of  in  modes  1  and  2  indicates  an  unstable 

node.  The  spatial  instability  modes  are  of  greatest  practical 
interest  and  most  of  the  eigenvalue  computations  described  in 
Section  3  are  using  these  modes. 

BMT  Ltd.  have  also  supplied  two  further  computer  programs 
known  as  BMTMF2  and  BMTST6  C63,  which  are  similar  to  NMIAL514  and 
NMIHE40S  but  deal  with  the  sixth-order  system  derived  by  Lowell  & 
Reshotko  C23,  which  Includes  fluctuations  in  all  the  flow 
variables.  They  found  that  buoyancy  effects  are  negligible.  The 
programs  only  deal  with  the  temporal  instability  case  and  a  few 
numerical  results  are  presented  in  Appendix  E. 


3.  NUMERICAL  RESULTS 

This  section  describes  test  computations  on  heated  and 
unheated  flat  plates  with  constant  wall  temperature.  The  cases 
chosen  were  those  used  by  Lowell  &  Reshotko  [23,  corresponding  to 
an  ambient  temperature  of  approximately  15.6°C  (60°F)  and  wall 
temperatures  of  approximately  15.6°C  (60#F),  32.2®C  (90°F), 

65.6*C  (150*F)  and  93.3*C  (200#F).  The  Reynolds  number  range  was 

from  400  to  20,000,  based  on  the  displacement  thickness  S*. 

This  is  typical  of  the  range  required  in  transition  calcul¬ 
ations,  since  S  is  generally  much  smaller  than  other  charac¬ 
teristic  length  scales,  such  as  body  length  or  diameter. 

Full  details  of  the  mean  flow  computations  are  presented  in 
Appendix  C.  The  mean  flow  profiles  were  converted  to  profiles  of 
dimensionless  velocity  and  its  second  derivative  (expressed  as 

functions  of  y/6*)  using  a  modified  version  of  program  NMI- 
AL514.  The  converted  profiles  were  input  to  a  modified  version  of 
program  NMIHE40S  to  calculate  the  eigenvalues. 

Lowell  &  Reshotko  C23  presented  the  results  of  their  eigen¬ 
value  calculations  in  terms  of  neutral  stability  curves,  which 

are  contour  plots  of  curves  of  ■  0  as  a  function  of  Reynolds 

number  R  and  «r.  To  obtain  comparisons  with  their  work. 


similar  plots  were  obtained  using  a  procedure  described  in 
Appendix  D. 

(a)  Unheated  Plates 

Some  tabular  comparisons  of  unheated  flat  plate 
eigenvalues  are  described  in  Appendix  0.  These  show  the 
compound  matrix  method  to  be  both  well-conditioned  numeri¬ 
cally  (only  single  precision  using  32-blt  arithmetic  is 
required)  and  capable  of  giving  eigenvalues  to  high  accu¬ 
racy.  Figure  1  shows  the  effect  of  the  number  of  inte¬ 
gration  steps  on  the  neutral  stability  curve.  The  results 
show  that  800  integration  steps  are  required  to  obtain 

accurate  eigenvalues  over  the  whole  range  of  &r  and  R. 

The  results  for  800  steps  agree  with  those  of  Lowell  & 
Reshotko.  The  eigenvalue  calculation  at  high  values  of 

oir  R  appears  to  be  sensitive  to  the  number  of  steps. 

At  these  values,  it  is  also  particularly  sensitive  to 

the  initial  guesses  for  and  w  (see  Appendix  D) . 

An  alternative  eigenvalue  solver  for  the  unheated  case, 
supplied  by  ARE  (Portland),  was  also  tested.  This  was 
originally  written  by  NPL  C7],  and  uses  a  matrix  method  with 
finite-difference  approximations  to  the  derivatives  of  the 
Orr-Sommerfeld  equation.  The  results  using  this  method  are 
presented  as  Figure  2.  There  was  less  variation  in  the 
solutions  obtained  using  100,  200  and  400  steps  than  with 
the  compound  matrix  method,  but  at  100  and  200  steps  the 
minimum  critical  Reynolds  numbers  were  calculated  to  be  557 
and  528  as  opposed  to  the  accepted  value  of  519  C83.  The 
solution  for  400  steps  using  this  method  and  all  the 
compound  matrix  solutions  for  200  steps  and  above  gave 
answers  within  1%  of  the  accepted  value.  However,  the 
accuracy  of  the  NPL  method  may  still  be  satisfactory  in 
engineering  calculations  of  flow  transition.  At  large 
numbers  of  steps  (greater  than  about  400),  the  method 
appears  to  suffer  from  ill-conditioning.  It  is  implemented 
in  double  precision  (64-blt  arithmetric) . 

The  effect  of  varying  the  freestream  boundary  is  shown 
in  Figure  3.  Using  800  integration  steps,  values  for 

of  10.0,  7.5  and  5.0  were  tried.  The  results  for  10.0  and 
7.5  were  almost  identical,  and  those  for  5.0  were  slightly 
different. 

(b)  ttefttg4.JBlAt.eg 

Some  numerical  results  are  presented  in  tabular  form 
and  discussed  in  Appendix  E.  The  effect  of  the  number  of 
integration  steps  on  the  neutral  stability  curve  was 
examined  for  the  cases  of  wall  temperature  32.2*C  (90*F)  and 
93.3*C  (200*F)  with  an  ambient  temperature  of  15.6*C  (60°F). 
The  freestream  boundary  in  the  mean  flow  calculation  was  at 
r\  ■  10  and  the  fluid  property  relations  were  those  of 
Lowell  &  Reshotko  C23.  In  the  former  case  (Figure  4), 
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FIG.  1  EFFECT  OF  NUMBER  OF  INTEGRATION  STEPS  ON 
NEUTRAL  STABILITY  CURVE  CALCULATED  USING 
COMPOUNO  MATRIX  METHOD  ( UNHEATED  PLATE) 
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FIG  2  EFFECT  OF  NUMBER  OF  INTEGRATION  STEPS  ON 
NEUTRAL  STABILITY  CURVE  CALCULATED  USING 
NPL  METHOD  171  (UNHEATED  PLATE) 


ARE  TM  (UHA)MIU 


-»  - 


MAL  PART  OF  WAVE  NUMBER 


800  integration  steps  are  sufficient  for  practical  calcu¬ 
lation  purposes,  but  in  the  latter  case  (Figure  5),  at  least 
1600  steps  are  required.  Figures  4  and  5  also  show  the 
nunerical  results  of  Lowell  &  Reshotko.  who  used  an  ortho- 
normalization  method.  At  the  lower  temperature  difference, 
their  results  agree  well  with  those  using  the  compound 
matrix  method  with  800  integration  steps.  At  the  higher 
temperature  difference,  however,  their  results  lie  somewhere 
between  those  using  800  and  1600  steps.  According  to 
Reference  2,  they  used  1000  steps,  which  the  present  results 
suggest  may  not  be  enough.  Their  values  of  dimensionless 
wall  shear-stress  and  temperature  gradient  agree  with  those 
in  the  present  investigation  (see  Appendix  C).  The  effect  of 
changing  the  freestream  boundary  is  shown  in  Figure  6 
corresponding  to  the  largest  overheat;  the  results  are 
similar  to  those  for  the  unheated  case  mentioned  above. 

The  effects  of  using  different  fluid  property  relations 
are  shown  in  Figures  7  and  8,  corresponding  to  small  and 
large  overheats  respectively.  Two  sets  of  relations  were 
tried,  those  used  by  Lowell  &  Reshotko  C23  and  those  of 
Kaups  &  Smith  C33  (for  details,  see  Appendix  B).  The 
assumption  of  constant  density  was  tried  in  the  Lowell  & 
Reshotko  formulae  and  made  very  little  difference  to  the 
calculated  eigenvalues. 

Figure  9  shows  the  effect  of  different  wall  tempera¬ 
tures  on  the  neutral  stability  curve.  These  results  were 
obtained  using  the  Kaups  &  Smith  fluid  property  relations 
(for  comparison  with  those  in  Reference  2),  with  variable 
density.  The  mean  flows  were  computed  using  an  accurate 
fourth-order  method  (Method  2,  Appendix  C).  Similar  plots 
were  obtained  using  different  fluid  property  relations 
(including  one  case  with  constant  density)  and  the  results 
(when  compared  graphically)  were  almost  indistinguishable 
from  those  in  Figure  9.  It  should  be  noted  that  the 
Reynolds  number  is  on  a  logarithmic  scale.  Further  results 
where  the  mean  flows  were  computed  using  a  second-order 
finite-difference  method  (Method  4,  Appendix  C)  were  also 
very  similar.  The  results  obtained  by  Lowell  &  Reshotko  are 
also  shown.  A  plot  of  minimum  critical  Reynolds  number 
against  overheat  (T  -  Tw)  is  presented  as  Figure  10.  This 

shows  the  data  collected  by  Lowell  &  Reshotko,  to  which  the 
present  numerical  results  have  been  added.  The  optimum 
overheat  is  seen  to  be  about  40°C. 

A  few  numerical  results  were  obtained  using  the  sixth- 
order  formulation  C63.  These  are  limited  to  the  temporal 
instability  case  (Mode  0)  with  a  wall  temperature  of  32.2°C 
(90*F)  and  are  described  in  Appendix  E.  Comparison  with  the 
fourth-order  formulation  shows  differences  in  eigenvalues 
which  are  nontrivial  but  small.  The  resulting  neutral 
stability  curves  (Figure  11)  correspond  with  those  of  Lowell 
&  Reshotko  C23.  It  should  be  noted  that,  on  the  neutral 
stability  curves  only,  the  spatial  and  temporal  eigenvalues 
coincide . 
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FIG  7  EFFECT  OF  DIFFERENT  FLUID  PROPERTIES  ON 
NEUTRAL  STABILITY  CURVE  OF  HEATED  PLATE 
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IFFERENT  FLUlO  PROPERTIES  ON  NEUTRAL  STABILITY  CURVE 


MESENT  RESULTS 
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FIG.  10  EFFECT  OF  OVERHEAT  ON  MINIMUM  CRITICAL  REYNOLDS  NUMBER 
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FIG.  11  COMPARISON  OF  NEUTRAL  STABILITY  CURVES 
USING  FOURTH -AND  SIXTH -ORDER  APPROACHES 
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*.  CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER  RESEARCH 

The  coapound  matrix  Method  gives  eigenvalues  which  generally 
agree  with  previously  published  results.  Different  expressions 
for  the  variation  of  fluid  properties  with  teaperature  give 
nontrivial  but  small  differences  In  calculated  eigenvalues,  as 
does  the  calculation  of  the  aean  flow  using  second-  and  fourth- 
order  numerical  methods.  Calculations  for  heated  flat  plates 
with  constant  wall  temperature  suggest  (In  accordance  with 
previous  work)  that  the  optimum  overheat  Is  about  40*C  and  that 
further  heating  has  a  destabilising  effect  on  the  boundary  layer. 

The  sixth-order  formulation,  where  fluctuations  of  all  the 
flow  variables  are  included,  gives  eigenvalues  which  are  only 
slightly  different  from  those  obtained  using  the  fourth-order 
formulation  based  on  the  modified  Orr-Sommerfeld  equation.  It  is 
not  practicable  to  incorporate  the  method  in  a  transition 
calculation  procedure  based  on  spatial  instability,  but  further 
eigenvalue  calculations,  based  on  temporal  instability,  could  be 
undertaken  if  thought  desirable. 

The  compound  matrix  method  requires  a  large  number  of 
integration  steps  (typically  BOO  to  1600)  to  obtain  accurate 
eigenvalues.  BMT  Ltd  C83  have  proposed  an  alternative  eigenvalue 
finder,  based  on  a  panel  method,  which  requires  relatively  few 
integration  steps.  They  demonstrated  that  the  method  works  using 
a  Blaslus  (unheated  flat-plate>  velocity  profile  in  the  temporal 
instability  case.  Further  research  would  be  necessary  for  a 
spatial  stability  implementation  to  be  used  in  a  transition 
prediction  method  and  this  might  be  useful  if  a  large  number  of 
heated  body  calculations  were  required. 

The  compound  matrix  method  is  suitable  for  inclusion  in  a 
method  for  predicting  boundary-layer  transition  on  axlsymmetric 
heated  bodies.  The  fluid  property  relations  for  fresh  water 
proposed  by  Lowell  &  Reshotko  C2J  are  suitable  for  use  in  such  a 
method.  The  assumption  of  constant  density  does  not  introduce 
any  appreciable  errors  and  it  can  be  used  to  simplify  an  axi- 
symmetrlc  calculation  procedure  and  to  reduce  the  computational 
effort.  A  second-order  finite-difference  boundary- layer  calcu¬ 
lation  method  could  be  used  to  determine  the  mean-flow  profiles. 

Although  it  simplifies  the  problem,  the  assumption  of 
constant  surface  temperature  is  physically  unrealistic.  To 
maintain  a  constant  surface  temperature  on  a  flat  plate  requires 
a  non-uniform  distribution  of  heat  flux  which  makes  the  mean  flow 
two-dimensional  rather  than  one-dimensional.  It  is  planned  for 
future  computations  to  prescribe  either  the  surface  teaperature 
or  the  surface  heat-flux  which  is  allowed  to  vary  with  streamwlse 
position. 
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APPENDIX  A 


Derivation  of  the  mean  flow  equations  for  a  heated  body 


This  Appendix  presents  the  equations  of  motion  of  a  heated 
two-dimensional  body  and  a  heated  flat  plate.  Let  x  represent 
the  co-ordinate  alonq  the  body  surface  and  y  the  co-ordinate 
normal  to  the  surface,  with  u(x,y)  and  v(x,y>  the  respective 
velocity  components  and  T(x,y)  the  absolute  temperature.  The 
fluid  properties  (density  p,  dynamic  viscosity  p,  specific  heat 
Cp  and  thermal  conductivity  k)  are  assumed  to  be  functions  of 

temperature  only.  The  usual  boundary-layer  assumptions  are  made, 
namely,  that  streamwlse  derivatives  are  much  less  than  cross- 
stream  ones.  The  equations  of  motion  are  C33:- 


Continuity: 

Ijlpul 

+  fy(PV>  *  0 

(Al) 

Momentum : 

(  3u 

pr  c 

* 

dU 

■  P-  ue  + 

iy(M  1^} 

(A2> 

Energy: 

pCp 

{“H  * 

V  m  =  i-  Ik 

ay/  ay  r 

aT\ 

ay/ 

(A3) 

In  the  momentum  equation  (A2>,  the  buoyancy  term, 

g*  <P  ~  P„> 

where  gx  is  the  acceleration  due  to  gravity  in  the  freestream 

direction,  has  been  neglected.  This  is  equivalent  to  assuming 
that  the  flow  Is  horizontal.  In  the  energy  equation  (A3),  the 
compression  work  term 


has  also  been  neglected,  since  (-1/p) (dp/dT)  (the  coefficient 
of  thermal  expansion  of  water)  Is  very  small  and  the  whole  term 
Is  also  small.  The  viscous  dissipation  term 

p  Ou/3y)2 

in  the  energy  equation  has  been  neglected,  since  Schlichting  C103 
has  shown  that  It  is  negligible  unless  the  Eckert  number 


E  ■  u2  /  ( C  T  ) 

m  p  “ 

Is  of  the  order  of  unity.  Ihe  value  of  E  In  the  present  appli¬ 
cation  Is  of  the  order  of  10~*. 
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The  velocity  boundary  conditions  are:- 


u  ■  0,  v  *  0  on  y=0 


(A4) 


u  -*  u@  as  y  -*  • 

Hie  teaperature  boundary  condition  at  the  wall  is  one  of  two 
options,  that  of  prescribed  wall  teaperature 

T  ■  Tw<x)  on  y  *  0  (A5) 


or  prescribed  wall  heat  flux 


-k  3T/8y  ■  4W<X>  on  y  *  0 


( A6) 


The  freestreaa  teaperature  boundary  condition  is 

T/T„  -*  1  as  y  -»  »  (A7) 


The  above  equations  (Al)  to  (A3)  with  the  boundary  con¬ 
ditions  (A4)  to  (A7)  are  transforaed  to  a  aore  convenient 
co-ordinate  systea  by  stretching  the  y  co-ordinate.  A  aodlf led 
Howarth-Dorotnitsyn  transformation  CIO]  is  used,  given  by 


n  * 


(A8) 

(A9) 


The  streaa  function  $  is  defined  in  the  usual  way,  as 


pu  =  Sdi/ 3y 
pv  *  -d$/9x 


(A10) 


and  a  dimensionless  streaa  function  f(n>  is  defined  by 


3f/3T)  *  u/ue 

(All) 

is  related  to  t  as 

follows : 

* 

»  <P„  ue  5)1/2  f 

(A12 ) 
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The  dimensionless  temperature  g(n>  Is  defined  by 


q  =  T/Tm 


(A13 ) 


The  equations  of  notion  are  transformed  to  C,  n  co-ordinates 
using  the  following  relations: 


where 


and 


L.  ,  i_ 

ax  ac 


*  "  ‘W  In 


,y 

I,  *  Op/Bx)  dy 

1  J0 

fy 

I-  =  p*  dy 

L  J0 


P 


a_ 

an 


(A14 ) 


(A15) 


The  continuity  equation  is  automatically  satisfied  by  the 
definition  of  the  stream  function.  After  transformation,  the 
momentum  and  energy  equations  become 

In lc  f">  *{r  - f'2} ' 5  (r  lr  - f''  If) 

{Pr-  Cp)  |^  «  5'l  *  f  »'  •  5  {''  |f  -  »'  H)  '«*> 


In  equations  (A16)  and  (A17),  dashes  denote  differentiation  with 
respect  to  n»  end 


C 


P  P  . 

P»  Mpp 


K 


£ 

P. 


k 


m 


L. 
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n  is  a  dimensionless  pressure  gradient  parameter  and  Pr^,  a 
freestream  Prandtl  number.  The  transformed  boundary  conditions 


are : 

f(0)  =  0, 

and  either  g(0) 

or  g'(0) 

As  f '  (r))  -»  1, 


f'(0)  =  0  (A18) 

»  gM  (A19) 

=  g^  ( A20 ) 

g ( rj)  ->  1  ( A21 ) 


Equations  (A16)  and  (A17)  along  with  the  boundary  conditions 
(A18 )  to  (A21)  describe  the  mean  flow  over  a  heated  two-dimen¬ 
sional  body.  In  the  case  of  a  flat  plate,  assuming  similarity  of 
the  velocity  and  temperature  profiles,  m  is  zero  and  (A16)  and 
(A17)  simplify  to 


|^(Cf">  +  iff"  =  0  ( A22 ) 

{Pr„  Cp}  1  (K  g')  +  |  f  g'  =  0  (A23) 

For  an  unheated  plate,  equation  (A22)  reduces  to  the  well- 
known  Blasius  equation 


f  "  ' 


+ 


f  " 


=  0 


(A24) 
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The  nean  flow  calculation  procedure  requires  values  for 
density,  viscosity  and  thermal  conductivity  and  their  first 
derivatives  and  values  for  specific  heat  capacity,  all  expressed 
as  functions  of  temperature.  The  numerical  procedure  used  by 
Lowell  &  Reshotko  C23  requires,  in  addition,  one  further  deriva¬ 
tive  of  each  quantity. 

Two  sets  of  fluid  property  relations  were  used  in  the 
present  investigations.  Those  proposed  and  used  by  Lowell  Si 
Reshotko  C23  are  as  follows: 

Specific  heat  capacity 

Cp/Cp0  =  2-13974  -  9-68137X10"3  T  +  2-68536xl0~5  T2 

-  2 -42139xl0-®  T3  (T  in  K)  (Bl) 


1_  d£  _  1  (T-a ) ( 3T-2b-a ) 

P0  dt  TT+dT  c 


+ 


+ 


exp(-f/T>  CT2  +  f ( T+d ) } 
Tz 


( B3 ) 


1_  d^£  _  2_  d£  1  |~(-6T+4a-2b) 

p0  dT2  "  P0  dT  (T+d)  L  c 


+  C ( f -2d )T+fd3  exp(-f/T)J  (B4) 


In  equations  (B2)  to  (B4),  T  is  in  °C  and 
a  =  3-9863,  b  =  288-9414,  c  =  508929-2,  d  =  60-12963 
e  =  0-011445  and  f  =  374-3 


Dynamic  viscosity 


log10(a  n0/P>  -  b(T-20-)d^9(J---?->; 


(B5) 


1  dyi  _ 
Ma  « 


r- 


+  2c(T-20)  -  log10(a  M0/p)> 

- eT3+!ri - 


<B6) 


L-  =  -  2  /E_£  +  L_  du\  +  _L 
m0  a?  W  aT/  mi 


(B7> 


In  equations  (B5)  to  (B7),  T  is  in  °C  and 

a  =  1-002,  b  =  1-37023,  c  =  8-36xl0“4,  d  =  109-0,  e  =  0-43429448 
Thermal  conductivity 


k/k0  =  -9-90109  +  0-1001982  T  -  1 -873892xl0-4  T2 

+  1-039570X10"7  T3 


(B8) 


Expressions  for  the  derivatives  of  Cp/CpQ  and  k/k0  can  easily 

be  obtained  from  equations  (Bl>  and  (B8).  In  the  above  ex¬ 
pressions. 


Cp0  =  4186-8  J  kg-1K-1  ,  p0  =  1000  kg  m  3, 
M0  *  0-001  N  s  nf2,  kQ  =  0-1  W  m_1K-1 


( B9 ) 


The  value  of  Cp0  differs  slightly  from  that  used  by  Lowell  & 

Reshotko,  but  this  has  been  shown  to  have  a  negligible  effect  on 
the  numerical  results. 


A  different  set  of  expressions  were  proposed  by  Kaups  & 

Smith  C33,  who  fitted  least-square  polynomials  to  the  normalised 
experimental  data  to  give 

C  2 

£ —  =  1-4833689  -  0-8072501  —  +  0-3289602  }  (B10) 

up,ref  Aref  l^refj 


where  a  =  35-15539,  b  =  -106-9718715,  c  =  107-7720376, 
d  =  -40-5953074.  and  e  =  5.6391948. 


Tref 

Mref 


•te)  [*  *  -fell 

-1-940589  +  5-2220185  j  -  2-693322  | 


( B13 ) 


+  0-4176176 


( B14 ) 


The  reference  values  in  the  above  expressions  are: 

Tref  =  273 ’ 15  K*  cp,ref  =  4218  J  k9~lK_1»  Pref  =  999-84  *9  m”3, 
Mref  =  0,0017936  N  8  n~2  and  kref  *  °-5537  H  m_1K_1  (B15) 


Figures  12  (a)  to  (d)  show  comparisons  of  the  two  sets  of 
formulae  for  Cp,  p,  p  and  k  compared  with  known  values 

C11.123  in  the  range  0-100°C  <273-15-373-15  K>.  There  is  little 
difference  between  them  for  p  and  k,  but  Lowell's  expressions 
agree  more  favourably  with  the  known  values  for  Cp  and  p.  The 

Kaups  &  Smith  formulae  were  designed  to  fit  the  data  over  a  much 
wider  temperature  range  than  that  of  Interest  in  the  present 
investigation.  The  formulae  of  Lowell  &  Reshotko  are  recommended 
for  use  in  future  investigations. 
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DENSITY  j  p  /  kg  SPECIFIC  HEAT  CAPACITY,  Cp/Jkf 


RG.  12(0)  COMPARISON  OP  DIFFERENT  FORMULAE  FOR 
SPECIFIC  HEAT  CAPACITY  WITH  DATA  FROM  TABLES 


FIG.  12(b)  COMPARISON  OF  DIFFERENT  FORMULAE  FOR 
DENSITY  WITH  OATA  FROM  TABLES 
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AMT*  (UHA)MtK 


THERMAL  CONDUCTIVITY,  k  /Wi  1  K  1  DYNAMIC  VISCOSITY,  IO*^/N*m 


FIG.  12(c)  COMPARISON  OF  DIFFERENT  FORMULAE  FOR 
DYNAMIC  VISCOSITY  WITH  DATA  FROM  TABLES 


FIG.  12(d)  COMPARISON  OF  DIFFERENT  FORMULAE  FOR 
THERMAL  CONDUCTIVITY  WITH  DATA  FROM  TABLES 
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ADC  TM(  UNAIHSU 


The  seen  flow  problem,  defined  by  equations  (A22)  and  (A23) 
with  boundary  conditions  (A18),  (A19)  and  (A21)  Is  a  two-point 
non-linear  boundary-value  problem.  Two  NAG  Library  C133  routines 
exist  for  this  type  of  problem,  one  called  D02GAF  (Method  1), 
based  on  a  finite-difference  type  method,  and  the  other  called 
D02HAF  (Method  2),  based  on  a  shooting  and  matching  method.  To 
use  these  methods,  equations  (A22)  and  (A23)  have  to  be  recast  In 
a  slightly  different  form,  viz: 


. . ft#  * 

-  TTX”” 

(Cl) 

■  -{£!?* 

( C2 ) 

Rather  than  calculating  out  to  infinity,  the  freestream  boundary 
conditions  are  applied  at  some  value  ne  which  is  specified  as 

a  numerical  parameter.  Most  of  the  present  calculations  were 
carried  out  using  *  10.  In  practice,  a  smaller  value  may 

be  acceptable  for  engineering  calculations.  Initial  estimates 
have  to  be  supplied  for  all  the  boundary  values,  i.e.  f.  f',  f''. 
g  and  g'  at  n  *  0  and  h  *  n  The  boundary  values  supplied 
were: 

f''(0)  *  1.0  for  unheated  boundary  layers , 
f ' ' ( 0 )  *  0.5  for  heated  boundary  layers, 
g'<0>  -  -0.05 

f(V  *  ’le  '  2 

f  "  < he )  »  0.0 
<? '  <  he )  «  0.0 

For  a  time,  the  author  did  not  have  access  to  the  NAG  Library, 
and  so  a  third  method  was  tried.  This  method  Is  the  one  used  by 
Lowell  &  Reshotko  C23  and  is  based  on  a  special  technique  used 
for  handling  asymptotic  boundary  conditions  described  by  Nacht- 
steln  &  Swlgert  C143.  Both  references  give  FORTRAN  llstlnas. 

This  method  only  requires  initial  estimates  for  f''(0>  and  g'(0) 
annd  the  values  used  for  these  were  the  same  as  those  used  In 
Methods  1  and  2. 

All  three  methods  were  used  to  obtain  mean  flow  profiles  for 
the  unheated  case  and  for  the  heated  cases  examined  by  Lowell. 

The  main  numerical  parameters  are  presented  in  Table  1.  It 
proved  difficult  to  obtain  solutions  using  Method  1,  with  ne 


having  a  value  of  10  In  all  the  heated  case*,  but  bom  solutions 
were  obtained  with  saaller  values  of  ne  giving  slightly 

different  answers.  The  values  of  f'*<0)  and  g'(0)  using  Methods 
2  and  3  agree  with  each  other  to  eight  declaal  places  and  they 
agree  with  the  values  of  Lowell  to  at  least  four  declaal  places. 
Method  2  was  found  easiest  to  use  as  the  output  profiles  could  be 
expressed  on  an  arbitrary  grid  of  points,  and  this  was  the  aethod 
used  for  aost  of  the  coaputations. 

The  saae  profiles  were  coaputed  using  a  second-order 
finite-difference  aethod  which  is  well  docuaented  in  the  litera¬ 
ture  C15.163  and  referred  to  as  Method  4  here.  The  wall-shear 
paraaeters  were  slightly  different  froa  those  using  fourth-order 
aethods  where  these  have  converged.  The  effect  of  these  differ¬ 
ences  is  shown  in  Appendices  0  and  E. 

Calculations  using  the  NAG  library  routines  were  carried  out 
in  double  precision  arlthaetic  (64-blt  word  length).  Most  of  the 
other  calculations  were  in  single  precision  arlthaetic  (32-bit 
word  length).  In  general,  double  precision  aade  no  appreciable 
difference  coapared  with  single  precision. 


Freestreas 


Aabient 

Mall 

Method 

Boundary 

f  "  (0) 

g’  ( 0) 

15.5556 

15.5556 

1 

10 

0.33205734 

2 

10 

0.33205734 

3 

10 

0.33205734 

Lowell 

10 

0.33205733 

4 

10 

0.33203962 

15.5556 

32.2222 

1 

10 

0,46640722 

-0.04118285 

2 

10 

0.46640722 

-0.04118285 

3 

10 

0.46640722 

-0.04118285 

Lowell 

10 

0.46639971 

-0.04117200 

4 

10 

0.46631959 

-0.04124339 

15.5556 

65.5556 

1 

5 

0.74604546 

-0.13595093 

2 

10 

0.74322894 

-0.13577964 

3 

10 

0.74322893 

-0.13577963 

Lowel  1 

10 

0.74319378 

-0.13574200 

4 

10 

0.74280715 

-0.13596509 

15.5556 

93.3333 

1 

8 

0.96324417 

-0.22613273 

2 

10 

0.96324381 

-0.22613270 

3 

10 

0.96324381 

-0.22613270 

Lowell 

10 

0.96317694 

-0.22606900 

4 

10 

0.96232432 

-0.22641809 

Table  I 


•CLT-rnr 

■rmr 


1 1 
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The  effect  of  different  fluid  property  relations  on  the  aean 

flow  parameters  and  the  displacement  thickness  S*  Is  shown  In 
Table  II.  All  computations  were  using  Method  2  with  an  ambient 
temperature  of  15.5556*C  and  the  freestream  boundary  at  n  ■  10. 
The  results  show  that  different  property  relations  give  small 
differences  in  the  numerical  parameters,  but  the  assumption  of 
constant  density  has  hardly  any  effect. 


Util  Fluid  Util  ifatAT  Hfcll  inett  nut  Displacement 

Temp  <»C)  ErgBCCtigl  flO)  g(0)  thickness  6* 


15.5556 

0.332057 

1.720768 

32.2222 

Lowell 

0.466407 

-0.041163 

1.514571 

Lowell  (const  p) 

0.466429 

-0.041212 

1.514638 

Kaups  &  Smith 

0.467292 

-0.041341 

1.514442 

65.5556 

Lowell 

0.743329 

-0.135780 

1.211281 

Lowell  (const  p> 

0.743271 

-0.135667 

1.211291 

Kaups  &  Smith 

0.743302 

-0.135909 

1.211802 

93.3333 

Lowell 

0.963244 

-0.226133 

1.037017 

Lowell  (const  p) 

0.963143 

-0.225632 

1.036866 

Kaups  &  Smith 

0.959619 

-0.225790 

1.038404 

r 


Detailed  nuwrlctl  results  for  unheated  elates 


RuMrictl  results  for  the  stability  of  the  Blaslus  boundary 
layer  corresponding  to  an  unheated  plate  have  been  published  In 
the  literature  by  Jordlnson  C 1 7]  for  the  case  of  spatial  lnsta- 
bllly  and  by  Davey  (referenced  by  Caster  C183)  for  the  case  of 
temporal  Instability.  For  comparison  with  these  results, 
eigenvalues  were  calculated  using  800  integration  steps  in  the 
former  case  and  3200  steps  in  the  latter  case,  with  the  free- 
streaa  boundary  at  n  *  10.  Table  III  shows  comparisons  for  the 
spatial  instability  case  and  Table  IV  for  the  temporal  insta¬ 
bility  case. 


Present  investigation  Jordlnson  U7J 


Remolds 

Freauencv 

Have  number  a 

Hftye_ 

number  a 

Number 

w 

Real 

lags 

Real 

Imag 

336 

0.1297 

0.308349 

0.007940 

0.30B4 

0.0079 

598 

0.1201 

0.307848 

-0.001897 

0.3079 

-0.0019 

998 

0.1122 

0.308596 

-0.005710 

0.3086 

-0.0057 

Table  III 

Comparlfqn 

of  eigenvalues  with 

tfagse  of 

Jordlnson  1 

Present  Investigation  Paver  [103 

Reynolds  Erequgpo  Have  number  a  Have  number  a 


Number 

w 

Real 

Imag 

Real 

Im»g 

500 

0.3 

0.11930379 

-0.00027998 

0.11930376 

-0.00027998 

1500 

0.2 

0.06312283 

0.00315660 

0.06312291 

0.00315663 

3000 

0.15 

0.04021918 

0.00278070 

0.04021919 

0.00278079 

Engineering  calculations  of  transition  on  axlsyaaetrlc 
bodies  are  noraally  carried  out  using  a  second-order  accurate 
finite-difference  eethod  CIS. 163  as  opposed  to  the  fourth-order- 
accurate  methods  used  specifically  for  a  flat  plate  calculation. 
To  see  whether  this  aade  any  significant  difference  to  the 
eigenvalues,  some  further  tests  were  carried  out  using  the 
Blaslus  boundary  layer  calculated  in  three  different  ways.  These 
were  (a)  a  fourth-order  calculation  using  Method  I  (see  Appendix 
C) ,  (b)  a  second-order  finite  difference  calculation  using  a 
Modified  version  of  the  prograa  described  in  Reference  IS.  with 
the  second  derivative  of  velocity  f'''(n>  calculated  using  a 
second-order  finite-difference  approxlaatlon.  and  (c>  with  the 
velocity  f‘(n>  calculated  using  a  second-order  net hod  and 
f '  "  <H>  calculated  exactly  froa  the  differential  equation  (A24). 
The  results  using  all  three  aethods  are  tabulated,  along  with 
Jordlnson's  results  In  Table  V.  In  all  cases,  the  freestreaa 
boundary  was  at  n  *  10  and  S00  Integration  steps  were  used. 


Remolds 

Frequency 

Method 

Have  number  a 

Muaber 

w 

Real  Iaag 

336 

0.1297 

(a) 

0.30835 

0.00793 

(b) 

0.30828 

0.00798 

<c) 

0.30842 

0.00794 

Jordinson 

0 . 3084 

0.0079 

598 

0.1201 

(a) 

0.30785 

-0.00190 

<b> 

0.30779 

-0.00184 

(c) 

0.30792 

-0.00190 

Jordinson 

0.3079 

-0.0019 

998 

0.1122 

(a) 

0.30860 

-0.00571 

(b) 

0.30856 

-0.00565 

( c  > 

0.30866 

-0.00572 

Jordinson 

0.3086 

-0.0057 

Thble  V  C«parl»pn--0f  eigenvalues  using  different  aean  flows 
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Cow pound  aatricea  WPL  Method 


Reynolds 

Freq 

No  of 

Nii»h«r  a 

Have  nuaber  a 

Nuaber 

w 

SiSER 

Real 

Iaaq 

Real 

Iaaq 

400 

0.02 

3200 

0.068109 

0.026870 

Failed  to 

converge 

1600 

0.068109 

0.026870 

0.073973 

0.028631 

800 

0.068109 

0.026871 

0.069476 

0.027321 

400 

0.068111 

0.026872 

0.068446 

0.026987 

200 

0.068122 

0.026895 

0.068200 

0.026915 

100 

0.068357 

0.027173 

0.068162 

0.026943 

400 

0.15 

3200 

0.355209 

0.005268 

0.363537 

0.003345 

1600 

0.355209 

0.005268 

0.357208 

0.004773 

800 

0.355209 

0.005268 

0.355703 

0.005160 

400 

0.355222 

0.005265 

0.355331 

0.005306 

200 

0.355379 

0.005206 

0.355235 

0.005538 

100 

0.356780 

0.004176 

0.355187 

0.006365 

1000 

0.01 

3200 

0.279B28 

-0.007287 

0.290591 

-0.009000 

1600 

0.279828 

-0.007287 

0.282377 

-0.007763 

800 

0.279832 

-0.007289 

0.280454 

-0.007387 

400 

0.279884 

-0.007313 

0.279976 

-0.007224 

200 

0.280443 

-0.007734 

0.279832 

-0.006920 

100 

0.282826 

-0.01344B 

0.279681 

-0.005818 

2200 

0.15 

3200 

0.390733 

0.039061 

0.402942 

0.035210 

1600 

0.390736 

0.039060 

0.393718 

0.038130 

800 

0.390776 

0.039051 

0.391486 

0.038786 

400 

0.391265 

0.038833 

0 . 390966 

0.038816 

200 

0.394557 

0.034749 

0.390944 

0.038324 

100 

0.382908 

0.011654 

0.391034 

0.036591 

10000 

0.03 

3200 

0.135634 

-0.011900 

0.156262 

-0.011298 

1600 

0.135645 

-0.011905 

0.140224 

-0.012065 

800 

0.135774 

-0.011985 

0.136744 

-0.011925 

400 

0.136942 

-0.013373 

0.135896 

-0.011769 

200 

0.134532 

-0.027835 

0.135639 

-0.011321 

100 

0.101152 

-0.032060 

0.135369 

-0.009668 

10000 

0.15 

3200 

0.389878 

0.036309 

0.400500 

0.032976 

1600 

0.389923 

0.036298 

0.392460 

0.035479 

800 

0.390448 

0.036039 

0.390498 

0.036111 

400 

0.393632 

0.031597 

0.389957 

0.036295 

200 

0.379159 

0.010141 

0.389628 

0.036469 

100 

0.310466 

0.036128 

0.388832 

0.037578 

Table  VI  Effect  of  nuaber  of  Integration  steps  on  eigenvalues 
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Some  comparisons  were  made  of  the  number  of  integration 
steps  on  the  eigenvalues.  The  comparisons  were  made  using  both 
the  compound  matrix  method  and  a  much  older  finite-difference 
type  matrix  method  developed  at  NPL  C173,  which  was  supplied  by 
ARE  (Portland).  Eigenvalues  were  generated  starting  with  a 
Blaslus  profile  at  101  points  and  converting  the  data  to  the 
required  number  of  integration  steps  using  the  modified  version 
of  program  NMIAL514.  The  freestream  boundary  was  at  n  *  10. 

The  results  are  presented  in  Table  VI  above.  They  show  that  the 
compound  matrix  method  produces  eigenvalues  which  converge  to  a 
“true"  value  as  the  number  of  integration  steps  is  Increased.  At 
high  Reynolds  numbers,  800  steps  are  required  to  produce  reason¬ 
ably  accurate  eigenvalues.  The  NPL  method,  when  used  with  100  to 
800  steps,  gives  eigenvalues  which  do  not  change  very  much  with 
the  number  of  steps.  The  eigenvalues  generally  lie  within  1%  of 
the  "true"  value,  except  near  the  region  of  neutral  stability 

where  is  small.  Above  800  steps,  however,  the  method  appears 

to  suffer  from  ill-conditioning.  It  is  implemented  in  double 
precision  (64-bit  word  length). 

Some  numerical  problems  were  encountered  when  using  the 
compound  matrix  method.  At  high  Reynolds  numbers  and  high  wave 
numbers,  premature  convergence  occured,  owing  to  a  floating-point 
underflow  condition.  The  problem  could  be  cured  by  running  in 
double  precision  with  a  special  option  known  as  G_FLOATING,  where 
the  smallest  positive  floating-point  number  is  approximately 

10-308^  opposed  to  single  precision  and  conventional  double 

precision  where  the  corresponding  number  is  approximately  10-38. 
However,  this  required  twelve  times  the  cpu  time  compared  with 
the  single  precision  version.  The  errors  arising  from  this 
problem  were  generally  in  the  sixth  place  of  decimals,  and  only 

occurred  in  cases  when  was  relatively  large  and  positive, 

hence  being  of  no  real  significance  in  a  practical  transition 
calculation. 

In  some  cases,  the  eigenvalue  calculation  was  sensitive  to 
the  initial  guess.  This  problem  occurred  when  using  single 
precision  and  double  precision  arithmetic,  both  with  and  without 
the  G_FLOATING  option.  Some  numerical  experiments  were  conducted 
using  a  Blaslus  profile  calculated  at  81  points  with  the  free¬ 
stream  boundary  at  r)  =  8  and  interpolated  to  801  points  (800 
integration  steps).  Both  modes  1  and  2  (spatial  instability 

modes)  were  tried  with  a  range  of  initial  guesses  for  »r  and 

(mode  1)  and  and  utp  (mode  2).  The  results  are  shown  in 

Figures  13  (mode  1)  and  14  (mode  2).  They  show  that  at  high 
Reynolds  numbers,  the  convergence  of  the  compound  matrix  method 
is  sensitive  to  the  initial  guess.  In  particular,  when  using 
mode  1,  to  have  any  chance  of  convergence,  the  initial  guess  for 

»r  should  err  on  the  low  side,  and  with  mode  2,  the  initial  guess 
for  5r  should  err  on  the  high  side.  Over  the  range  of  Reynolds 


cttco  0 


FIG.  13  GRAPH  SHOWING  AREA  OF  CONVERGENCE  FOR  TYPICAL  EIGENVALUE 
CALCULATIONS  USING  MODE  1 


y 


numbers  tested  (400  to  20000)  this  was  only  found  to  be  a  problem 
with  unheated  boundary  layers  at  Reynolds  numbers  above  1000. 

Although  the  eigenvalue  calculation  method  was  initially 
supplied  as  an  interactive  program,  many  of  the  computations 
involved  generating  a  series  of  eigenvalues  over  a  range  of 
Reynolds  number  and  wave  number.  Initially,  the  results  of  one 
eigenvalue  calculation  were  used  as  the  initial  guess  for 
another,  but  this  method  often  failed  at  high  Reynolds  numbers. 

To  overcome  this  difficulty,  the  following  procedure  was  adopted. 
In  mode  2,  used  for  most  of  the  calculations,  for  the  lowest 

Reynolds  number  and  the  smallest  value  of  a  ,  user-specified 
values  of  and  iUr  were  used  as  initial  guesses.  The  result¬ 
ing  value  of  o£^  was  used  as  the  initial  guess  at  the  next  value 
of  <xr  and  the  initial  guess  for  iUr  was  based  on  linear  extrapo¬ 
lation  with  af.  At  the  second  lowest  Reynolds  number,  for  each 
ixr,  the  corresponding  values  of  oi1  and  at  the  first  Reynolds 
number  were  used  as  initial  guesses.  At  higher  Reynolds  numbers, 
for  a  particular  ot^,  the  initial  guesses  for  and  were 

calculated  using  linear  extrapolation  with  ln(R).  This  pro¬ 
cedure  worked  well  for  all  the  heated  cases,  but  no  scheme  could 
be  found  which  worked  successfully  for  all  the  unheated  cases 
tested. 


-43- 


* 


IV 

u. 

o 


<  o 

UJ 

in  ►*> 

<  < 

UJ 

—  i 

IV 

< 

u. 

O  tt: 

o 

in  li¬ 
ar 

3  a 
o  uj 
*-  CD 

Z  X 

o  3 

O  X 

u.  m 

o  o 

-j 

I-  o 

o  X 

-I  V 

cl  m 

.  a: 


^  z 
0-  < 


AACTM(UHA)|tB14 


1 

j 

APPENDIX  E 

Detailed  numerical  results  for  heated  Plates 

The  effect  of  using  different  numbers  of  integration  steps 
on  a  few  typical  eigenvalue  calculations  using  node  2  is  shown  in 
i  Table  VII  below.  The  results  show  that  a  large  number  of 

j  integration  steps  are  required  to  predict  the  eigenvalues 

|  accurately.  In  particular,  to  obtain  accurate  values  in  unstable 

regions,  where  the  imaginary  part  of  the  wave  number  is  negative, 

I  at  least  1600  steps  are  required. 

The  effect  of  different  mean  velocity  and  temperature 
profiles  on  eigenvalue  calculation  is  shown  in  Table  VIII  below. 
Different  profiles  were  obtained  by  using  different  expressions 
for  the  fluid  property  variation  with  temperature.  Three  sets  of 
expressions  were  used,  being  those  of  Lowell  6  Reshotko  C23  with 
both  variable  and  constant  density,  and  those  of  Kaups  &  Smith 
C33.  These  profiles  were  obtained  for  each  value  of  wall 
temperature,  using  a  fourth-order  "exact"  solution  (Method  2  in 
Appendix  C).  A  fourth  set  of  profiles  was  obtained  using  a 
second-order  "approximate"  f inite-diff erence  solution  (Method  4 
in  Appendix  C),  with  the  Lowell  &  Reshotko  expressions  for 
fluid-property  variation  with  variable  density.  The  results  in 
Table  VIII  indicate  that  for  a  particular  wall  temperature,  small 
changes  in  the  mean-flow  profiles  (which  are  quantified  in 
Appendix  C)  lead  to  relatively  small  changes  in  the  eigenvalues. 

.  In  particular,  the  assumption  of  constant  density  hardly  affects 

the  eigenvalue  calculation.  An  errors  obtained  using  this 
assumption  are  likely  to  be  smaller  than  errors  obtained  using 
too  few  integration  steps. 

A  typical  plot  of  contours  of  a ^  as  a  function  of  R  and 

aT  is  presented  as  Figure  15.  In  this  case,  the  wall  temperature 

is  32.2°C  (90#F)  and  the  Lowell  &  Reshotko  property  variation  is 
used  with  800  integration  steps.  The  results  using  1600  steps 

show  only  very  slight  differences  in  the  unstable  region  (« 
negative ) . 

Figure  15  was  obtained  by  computing  a  table  of  eigenvalues 
for  a  few  values  of  and  R  (normally  eight  and  ten  respect¬ 
ively)  using  mode  2.  To  produce  a  smooth  contour  plot  and  to 
save  on  computer  time,  a  much  larger  table  of  eigenvalues  was 

generated  for  51  uniformly-spaced  values  of  and  51  uniformly- 

|  spaced  (on  a  logarithmic  scale)  values  of  R,  using  cubic  spline 

interpolation  (NAG  Library  routine  E01ACF  C133).  Contour  plots 
were  drawn  using  the  larger  table  using  the  SIMPLEPLOT  Mark  I 
graphics  package.  All  of  the  neutral  stability  curves  presented 
above  were  produced  in  a  similar  way.  In  general,  the  contour 
plots  obtained  were  smooth,  but  some  of  those  for  the  unheated 
case  did  have  slight  ripples  which  were  smoothed  on  the  final 
■  versions  of  the  figures.  The  interpolation  program  was  also  used 

to  determine  the  minimum  critical  Reynolds  number  below  which  all 
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disturbances  are  stable.  The  values  for  various  cases  are  in 
Table  IX  below.  They  are  discussed  further  in  Section  3. 


Nall 

Reynolds 

Number  of 

Navenumber 

Freauencv 

Temperature 

Number 

Integration 

a 

w 

(•C) 

R 

steps 

Real 

In  a  a 

32.2222 

400 

3200 

0.05 

0.024346 

0.016689 

1600 

0.024346 

0.016689 

800 

0.024346 

0.016689 

400 

0.024347 

0.016689 

2000 

3200 

0.10 

0.012078 

0.025248 

1600 

0.012078 

0.025248 

800 

0.012081 

0.025246 

400 

0.012109 

0.025215 

20000 

3200 

0.15 

-0.002720 

0.026628 

1600 

-0.002769 

0.026617 

800 

-0.003512 

0.026507 

400 

-0.011928 

0.026425 

65.5556 

400 

3200 

0.05 

0.032698 

0.020750 

1600 

0.032698 

0.020750 

800 

0.032698 

0.020750 

400 

0.032705 

0.020750 

2000 

3200 

0.10 

0.020908 

0.022470 

1600 

0.020908 

0.022470 

800 

0.020917 

0.022465 

400 

0.020998 

0.022398 

20000 

3200 

0.15 

-0.002655 

0.022820 

1600 

-0.002744 

0.022794 

800 

-0.004294 

0.022565 

400 

-0.020114 

0.023288 

93.3333 

400 

3200 

0.05 

0.048665 

0.017552 

1600 

0.048665 

0.017552 

800 

0.048666 

0.017552 

400 

0.048691 

0.017553 

2000 

3200 

0.10 

0.021585 

0.020538 

1600 

0.021587 

0.020538 

800 

0.021613 

0.020525 

400 

0.021804 

0.020348 

20000 

3200 

0.15 

-0.004155 

0.021116 

1600 

-0.004436 

0.021047 

800 

-0.009386 

0.020606 

400 

-0.028066 

0.024785 

Table  VII  Effect  of  number  of  integration  steps  on 
heated  flat-plate  eigenvalues 
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Wall 

Fluid 

Revnolds 

Wave 

Freauencv 

Temperature  Properties/ 

Number 

Number,  a 

w 

(•C) 

Mean  flow 

R 

Real 

Imaq 

32.2222 

(S00  Int 
steps) 

Lowell /2nd  order 
Lowell /4th  order 
Lowell  (const  p) 
Kaups  &  Smith 

400 

0.05 

0.024367 

0.024346 

0.024339 

0.024297 

0.016687 

0.016689 

0.016688 

0.016698 

32.2222 

Lowell /2nd  order 
Lowell/4th  order 
Lowell  (const  p) 
Kaups  &  Smith 

2000 

0.10 

0.012070 

0.012081 

0.012086 

0.012122 

0.025247 

0.025246 

0.025248 

0.025255 

32.2222 

Lowell /2nd  order 
Lowell/4th  order 
Lowell  (const  p) 
Kaups  &  Smith 

20000 

0.15 

-0.003488 

-0.003511 

-0.003504 

-0.003488 

0.026511 

0.026507 

0.026506 

0.026500 

65.5556 

(1600  lnt 
steps) 

Lowell/2nd  order 
Lowell/4th  order 
Lowell  (const  p) 
Kaups  &  Smith 

400 

0.05 

0.032621 

0.032698 

0.032690 

0.032807 

0.020642 

0.020750 

0.020766 

0.020832 

65.5556 

Lowell /2nd  order 
Lowell /4th  order 
Lowell  (const  p) 
Kaups  &  Smith 

2000 

0.10 

0.020859 

0.020909 

0.020923 

0.020932 

0.022467 

0.022470 

0.022467 

0.022464 

65.5556 

Lowell /2nd  order 
Lowell/4th  order 
Lowell  (const  p) 
Kaups  &  Smith 

20000 

0.15 

-0.002712 

-0.002744 

-0.002750 

-0.002782 

0.022797 

0.022794 

0.022793 

0.022797 

93.3333 

(1600  int 
steps ) 

Lowell /2nd  order 
Lowell /4th  order 
Lowell  (const  p) 
Kaups  &  Smith 

400 

0.05 

0.04B365 

0.048665 

0.048747 

0.048943 

0.017528 

0.017552 

0.017540 

0.017527 

93.3333 

Lowell /2nd  order 
Lowell/4th  order 
Lowell  (const  p) 
Kaups  &  Smith 

2000 

0.10 

0.021542 

0.021587 

0.021574 

0.021520 

0.020533 

0.020538 

0.020532 

0.020537 

93.3333 

Lowell /2nd  order 
Lowell/4th  order 
Lowell  (const  p) 
Kaups  &  Smith 

20000 

0.15 

-0.002663 

-0.004436 

-0.004461 

-0.004547 

0.028493 

0.021047 

0.021047 

0.021072 

Table  VIII  The  effect  of  different  fluid  property  relations  and 
different  mean  flows  on  heated  flat-plate  eigenvalues 
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Wall 

Number  Of 

Fluid 

Mean-flow 

Minlauw  critical 

Temperature 

(•C) 

Intecrration 

Steps 

Properties 

Calculation 

Reynolds  number 

15.5556 

600 

2nd 

order 

521 

4th 

order 

519 

32.2222 

800 

Lowell 

2nd 

order 

5142 

Lowell 

4  th 

order 

5136 

Lowell  (const  p) 

5147 

Kaups  &  Smith 

5199 

65.5556 

1600 

Lowell 

2nd 

order 

11345 

Lowell 

4th 

order 

11297 

Lowell  (const  p) 

11289 

Kaups  &  Smith 

11227 

93.3333 

1600 

Lowell 

2nd 

order 

8340 

Lowell 

4th 

order 

8281 

Lowell  (const  p) 

8240 

Kaups  &  Smith 

8101 

Table  IX  Effect  of  different  fluid  properties  and  different 


wean  flows  on  minimum  critical  Reynolds  number 


A  few  numerical  results  were  obtained  for  the  teaporal 
stability  case  (node  0)  using  the  sixth-order  program  BMTST6  and 
these  are  coapared  with  the  corresponding  fourth-order  results  in 
Table  X.  The  wall  teaperature  was  32.2 °C  (90*F),  the  freestreaa 
boundary  at  h  =  10  and  the  Lowell  fluid-property  relations  were 
used  with  800  integration  steps.  In  general,  the  results  show  a 
non-trivial  but  saall  difference  between  the  fourth-order  and 
sixth-order  approaches.  Neutral  stability  curves  (which  are  the 
sane  for  teaporal  and  spatial  instability)  are  coapared  in  Figure 
15,  which  shows  that  the  differences  between  the  two  cases  is  the 
saae  as  that  obtained  by  Lowell. 


Reynolds 

Wave 

Fourth- 

order 

Sixth-order 

Number 

Number 

Frequenc 

y ,  w 

Frequency,  u> 

R 

a 

Real 

Inaa 

Real  Imaq 

5000 

0.12000 

0.026885 

-0.000573 

0.026820 

-0.000705 

6000 

0.16C00 

0.037071 

0.000190 

0.036080 

0.000138 

7000 

0.12233 

0.025799 

0.000080 

0.025788 

0.000003 

8950 

0.10668 

0.021016 

0.000074 

0.021014 

0.000002 

14250 

0.17315 

0.033227 

0.000081 

0.033332 

0.000000 

40  0 

0.05000 

0.015543 

-0.007361 

0.015300 

-0.006302 

2000 

0.10000 

0.025579 

-0.003712 

0.024762 

-0.004312 

10000 

0.15000 

0.030417 

0.000640 

0.030477 

0.000591 

20000 

0.15000 

0.026593 

0.000683 

0.026732 

0.000524 

Table  X 


